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Abstract 



' Many networks of physical and biological interest are characterized by a long-range coupling 

' mediated by a chemical which diffuses through a medium in which oscillators are embedded. We 



considered a one-dimensional model for this effect for which the diffusion is fast enough so as to 



i be implemented through a coupling whose intensity decays exponentially with the lattice distance. 



In particular, we analyzed the bursting synchronization of neurons described by two timescales 
(spiking and bursting activity), and coupled through such a long-range interaction network. One 
of the advantages of the model is that one can pass from a local (Laplacian) type of coupling to 
a global (all-to-all) one by varying a single parameter in the interaction term. We characterized 
bursting synchronization using an order parameter which undergoes a transition as the coupling 
parameters are changed through a critical value. We also investigated the role of an external 
time-periodic signal on the bursting synchronization properties of the network. We show potential 
applications in the control of pathological rhythms in biological neural networks. 
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I. INTRODUCTION 



In many problems of physical and biological interest we consider nonlinear oscillators 
whose interaction is mediated by a substance which is secreted by the cells and diffuses 
along the inter-cellular medium, being absorbed by the cells. Moreover, the rate of secretion 
depends on the cell dynamics, as well as the rate of absorption. In this way the dynamics of 
the cells are effectively coupled by the diffusing substance, leading to a non-local coupling 
type which depends on the details of the diffusion process as well as the dynamics of the 
oscillators themselves, leading to a complex system which displays a wealth of dynamical 
behaviors, like periodic and chaotic regimes, bifurcations, crises, destruction of tori, among 
other features. 

An outstanding biological example of interaction mediated by a chemical is the collective 
behavior of brain cells responsible by the circadian rhythm. The circadian rhythm is a daily 
periodicity (roughly a 24 hr cycle) of physiological, biochemical, and behavioral processes in 
living beings It is produced in mammals by specialized cells (circadian clocks) belonging 
to the suprachiasmatic nucleus (SCN) of the anterior hypothalamus The SCN consists of 



multiple, single-cell circadian clocks which, when synchronized, produce a coherent circadian 
output that regulate overt rhythms js-^. 

The circadian master clock of the SCN is entrained by the daily light-dark cycle, which 
acts via retina-to-SCN neural pathways jsj. Hence, to obtain a coordinated circadian rhythm, 
the master clock cell must be coupled to the other cells in the SCN so as to synchronize 
them to its own rhythm. The chemical coupling between circadian clocks can be described 
by means of a chemical (7-aminobutyric acid, or GABA, for short) which is both secreted 
and absorbed by clock cells immersed in some the intra-cellular medium. The coupling, in 
this case, is non-local in the sense that it takes into account cells which are not necessarily 
close to each other. An extreme situation belonging to this general category is the so-called 
global coupling, in which each cell interacts with the average concentration of the chemical 
in the inter-cellular medium due to all the other cells 

Another situation in which this kind of coupling is potentially important is chemotaxis, 
which is the influence of chemicals on the motion of somatic cells, bacteria, and other single 



or multiple-cell organisms jsj. These chemicals diffuse in the environment of the cells and 
the corresponding concentration gradients are responsible for the movement of the cells. For 



example, bacteria find glucose by swimming towards regions of higher concentration of such 



food molecules in the environment jo]. Alternatively, they can also flee from poison (such as 
phenol) according to their local concentration. Moreover, sperm can move towards the egg 
during fertilization, thanks to concentration gradients of aminoacids, in a ligand/receptor 
interaction. 

The chemotactic ability of slime molds like Dictyostelium is responsible for starving amoe- 
bae to gather to form multicellular bodies [lo| . During most of their lives, these slime molds 
are individual unicellular protists living in similar habitats and feeding on microorganisms. 
In the absence of food, however, they release signal molecules (DIF-1, short for Differentia- 
tion Inducing Factor) into their environment, so that they can find other amoeba and create 
swarms. In other words, when a chemical signal is secreted, they assemble into a cluster 



which acts as a single organism 



A model for nonlocal coupling mediated by a diffusing chemical was proposed by Ku- 
ramoto, in which the equations governing the time evolution of the oscillators can be coupled 
by using the concentration of a substance which diffuses through the medium in which the 



oscillators are embedded 12|. If the chemical diffuses in a timescale much faster than the 
oscillator period the coupling, although involves virtually all oscillators like in th e g lobal 
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case, depends on the distance between oscillators in an exponentially decaying way 
This approach has been used in studies of cell interaction |l5j and neural oscillators [16]. In 
the latter example the chemical secreted and absorbed by the neurons is a neurotransmitter 
which mediates the coupling among neurons. This coupling, on its way, depends also on 
the dynamics of the individual neurons, since it determines the propagation of electrical 
impulses along neuronal circuits in the brain. 

This neuronal activity (i.e., the evolution of the action potential) presents a fast timescale 
characterized by repetitive spiking and a slow timescale with bursting activity, where neuron 
activity alternates between a quiescent state and spiking trains There are many models 
for this spiking-bursting behavior, comprising differential equations like the Hindmarsch- 
and discrete-time maps, as the Rulkov map 



Rose model 



An assembly of coupled bursting neurons exhibit many self-organized phenomena. We 
are particularly interested in bursting synchronization, for which the neurons burst at ap- 
proximately the same time, even though their spiking behavior may be not synchronized 
itself jl^. Thanks to the slow timescale we may regard each bursting neuron as a phase 



oscillator, on defining a bursting phase and a corresponding frequency (its time rate) 22|. 

On the other hand, neurophysiologists argue that bursting synchronization plays a key 
role in some pathologies like Parkinson's disease, essential tremor, and epilepsies 23|. There- 
fore, the question of how to suppress this synchronization acquires a practical importance 
in terms of the control of undesirable neuronal rhythms. There have been proposed some 
alternatives to implement such a control through deep brain stimulation techniques. One 
of these is the addition of a time-periodic external signal of small amplitude and given fre- 



quency 
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24| . We have considered such a scheme in scale- free networks 26| and regular 



lattices with a power-law coupling [27]. Another deep brain stimulation technique was pro- 
posed by Rosenblum and Pikowsky, and consists of a time-delayed feedback signal applied to 
specific cortex areas 24, l28|]. We have recently studied this technique in scale- free networks 



29|. 



In this work we propose to study the control of bursting synchronization in a neuronal 
network with a long-range coupling mediated by a diffusing substance, using the Kuramoto 
model in the adiabatic limit of a rapidly diffusing chemical 12Nl4|. We shall use a time- 
periodic external signal and a time-delayed feedback signal so as to suppress bursting syn- 
chronization whenever it occurs, and study this phenomenon from the point of view of our 
coupling model. 

This paper is organized as follows: in Section 2 we outline a model for long-range coupling 
intermediated by a diffusing substance. In Section 3 we particularize the model to one- 
dimensional oscillator chains and coupled map lattices. In Section 4 we consider the bursting 
dynamics as described by the Rulkov map, and a long-range coupled network of Rulkov 
neurons. In Section 5 we investigate the bursting synchronization of coupled Rulkov neurons 
according to the nonlocal coupling model. Section 6 deals with the control of synchronized 
bursting rhythms by application of an external time-periodic driving and a delayed feedback 
signal. Our conclusions are left to the last Section. 



II. LONG-RANGE COUPLING MEDIATED BY A DIFFUSING SUBSTANCE 



Chemical coupling mediated by a diffusing substance can be mathematically described 
by a model proposed by Kuramoto leading to nonlocal couplings 12|. In this model the 
state variables of each oscillator influence the secretion of a chemical substance obeying a 
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diffusion equation [l3Hl6| . The rate of absorption depends on the local concentration of this 
substance at each cell position. In the following we will deal with two classes of vectors, 
which are represented with a different notation: (i) positions r in a (i-dimensional Euclidean 
space, to which the oscillators belong; (ii) state variables X = {xi,X2, ■ ■ ■Xm)'^ in a M- 
dimensional phase space of the dynamical variables characterizing the state of the system 
at a given time t. 

There are N oscillator cells located at discrete positions fj, where j = 1,2, - ■■ N, in the 
(i-dimensional Euclidean space; and Xj is the state variable for each oscillator, whose time 
evolution is governed by the vector field F(Xj) [Fig. [T]. The oscillators are not supposed to 
be identical, though, for they can have slightly different parameters. 

We suppose that the time evolution is affected by the local concentration of a chemical, 
denoted as A{f,t), through a time-dependent coupling function g: 

"^^^ F(X,) + g(A(f,t)), (1) 



dt 

whereas the chemical concentration satisfies a diffusion equation of the form 

^M(r) ^ ^ DV'A{f,t) + J2h{X,)6{f- n), (2) 

k=l 

where e <^ 1 is a small parameter representing the fact that diffusion occurs in a timescale 
faster than the intrinsic period of individual oscillators; rj is a phenomenological damping 
parameter, and D is a diffusion coefficient. The diffusion equation above has a source term 
h which depends on the oscillator state at the discrete spatial positions fj: this means that 
each oscillator secrets the chemical with a rate depending on the current value of its own 
state variable. 

According to Ref. 



12| we assume that the diffusion is very fast, compared with the 
oscillator period, such that we may set eA = 0. This makes the concentration of the 
mediating chemical to relax immediately to a stationary value that can be written in the 
following form: 

TV 

A{f,) = J2^i^J-^k)h{X,), (3) 

k=l 

where a{fj — fk) is a Green function (or an interaction kernel), which is the solution of 

{r] - DV^)a{fj -f) = 6{fj). (4) 



X 



FIG. 1. Schematic figure of spatially distributed phase oscillators. 



Since we have eliminated adiabatically the concentration of the diffusing chemical, on 
substituting ([3]) into Eq. ([2]) we can obtain an equation expressing the nonlocal coupling in 
the adiabatic approximation 

^ = F(X,) + - ^"'^)MX,) j . (5) 

If g is a linear function of Xj (but not necessarily of the positions fj) we can write 

^ = F(X,) + f;a(f, - f,)g(MX,)). (6) 

fc=i 

in such a way that we distinguish among some cases of interest: (i) the linear coupling, for 
which 

g(/i(Xfe)) = AXfc, (7) 

where A is a M x M matrix indicating which variables of the oscillators are coupled; (ii) 

the future coupling, where 

g(/i(X,)) = AF(Xfc), (8) 
and (iii) the nonlinear coupling, such that 

g(/i(Xfc)) = AH(Xfc), (9) 
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where H is a nonlinear function of its arguments. 

The equation for hnear couphng, for instance, is thus 

N 



dt 



F(X,) + 5^a(r-;-ffc)AX, 



(10) 



k=l 



As an example, let us consider each cell as undergoing a time evolution governed by the 
Hindmarch-Rose equations jisl, which also describe neurons with spiking and bursting ac- 
tivity. For this system we have M = 3 and 



X 



y 



^ y + ax^ — a;^ — z + / ^ 

1 — hx^ — y 
y r[s{x -x)- z] J 



(11) 



where a, b, I, r, s, and x ^i-re model parameters. 

Coupling these equations through the x-variable amounts to choose 



£00 






(12) 



where e is the coupling strength; such that the equation flTU]) for the linear coupling reads, 
in this case 

N 

- Zj + I + a{fj - fk)xk, 



k=l 



yj = l-hx]~yj, (j = l,2,...A) 

Zj = r[s{xj - x) - Zj]. 



(13) 



V 



If we were to couple these equations through the variable, we would have to use a different 
matrix, namely 

^0 0^ 

e , (14) 
^ 

and so on. The corresponding equation for future coupling is obtained simply by changing 
Xj by F(Xj) in 

As an example of nonlinear coupling, let us consider the case of M = 1, for which Xj is 
a single phase 9j G [0, 2tt), the vector function F(Xj) being the corresponding frequency coj 



(different for each oscillator, in general). In this case A reduces to a scalar coupling strength 
e and we can choose as the nonlinear coupling function 



H(Xfc) = sin(e,- - Ou 



(15) 



yielding a nonlinearly coupled Kuramoto model 30|, |31| 



N 

9j = Uj + e (T{fj — Vk) sin( 

A;=l 



Ok). 



(16) 



From the Fourier transform of Eq. (jlj), the interaction kernel can be written as 

aiVj - r) = ^ / d q — ^. 



(17) 



If the system is isotropic, the kernel becomes a function of the distance R = |rj — r| only, 
and can be expressed as 

^Ce-^^, iid=l, 
CKoi^R), ifc/ = 2, (18) 



a{R) = < 



C- 



7R 



if d = 3 



where Kq is the modified Bessel function of the second kind and order 0, the constant 7 is 
the inverse of the coupling length and is given by 
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(19) 



and the constant C is determined from the normalization condition 



d'^ra(r) = 1. 



(20) 



III. ONE-DIMENSIONAL LATTICES WITH LONG-RANGE COUPLING 



The most elementary application of the model of nonlocal coupling mediated by a diffusing 
chemical is a one-dimensional regular lattice of (an odd number) fixed and equally spaced 
oscillators with non-local interactions given by the kernel (!T8|) for d = 1. Assuming that the 
distance between consecutive lattice sites is a constant A, and supposing periodic boundary 
conditions 



±N' 



N' 



N -1 



(21) 



we can write \rj — fk\ = {j — k)A = iA by changing the summation index from k = 1,2, . . . N 
to i = j — k, such that i = ±1, ±2, . . . ± N'. This leaves us with 2N' + 1 = sites, each of 
them with a distance £ from any site j. The Green function is thus 

a{\r^-r\) = Ce-"'^'. (22) 

On excluding self-interactions, or the coupling of any site with itself, we replace the sum 
over k in Eq. (fTOj) with two sums, one over £ = 1,2, . . . N' and other over £ = —1, —2, . . . — N'. 
In the latter sum we can change index again m = —£ and then replace mhj £ since they are 
dummy indexes. Hence the first sum considers k = j — £ and the second one k = j + £, such 
that we can group them together into a single summation, giving for the linear coupling the 
following expression 

N' 

^ = F(X,) +CJ2 e-^^'A (X,_, + X,+,) . (23) 

i=i 

The normalization constant is determined from (|20|) which, in this case, becomes a sum- 
mation rather than an integral 

N 

^a(|r,-r|) = l. (24) 

k=l 

Making the same changes of index as in the previous paragraph we obtain 

-1 



C 



N' 

7A£ 



e 
1 



(25) 



On returning to our previous example of linearly coupled Hindmarch-Rose equations the 
x-coupled system is then 

N' 



Xj = Uj + ax^j — — Zj + I + eC e "'^^{xj-i + Xj 

1=1 

y, = l-hx]~y„ {3 = 1,2,. ..N) (26) 

Zj = r[s{xj-x) - Zj]. 

whereas, in the one-dimensional case, the nonlinearly coupled Kuramoto model fll6l) becomes 

N' 

9j = Uj + eCj2 e"^"^^ [sin(% - ^j-^) - sin(% - 9j+e)] . (27) 



e=i 



It is instructive to explore the limiting cases of the nonlocal coupling. If 7 goes to zero 
then, 

1 1 

^=2N' = iV^' ^^^^ 
and we have a global type of coupling 

^ = F(X,) + AX, (29) 

where each oscillator interacts with the mean field of all sites (except itself) 

N 

I 

X 



^— E (30) 



N _ 

In the limit of 7 large, the exponentials in the Green function decay very fast with 
the lattice distance such that only the term with £ = 1 contributes significantly to the 
summations. This gives for the normalization constant 

C « ^ (31) 

and the coupling term takes into account effectively only the nearest neighbors of a given 
site 

^ = F(X,) + ^A(X,_i + X,^0 (32) 
which is the usual local (or laplacian) coupling. 

IV. NONLOCALLY COUPLED RULKOV NETWORKS 

In this and the forthcoming sections we present, as an application of nonlocally couplings 
mediated by a diffusing substance, a model of coupled bursting neurons. Depending on 
the phenomena we are interested to investigate, the mathematical description of biological 
neurons may require models ranging from dozens of complicated differential equations to 



simple integrate-and-fire one-dimensional models j32|. In the following we are going to study 
collective phenomena involving the spiking and bursting activity of neurons, for which there 
is a slow bursting modulating the fast action-potential spiking. 

A. Local dynamics 



Since we are interested on finding the conditions for an assembly of neurons to burst 
synchronously, in particular the coupling parameters necessary to achieve this goal, the fine 
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details of the neuron dynamics are not essential to the network model. Hence, instead of 
continuous-time models like the Hindmarch-Rose model, which can be slower to simulate 
computationally, we choose to use discrete-time models instead, which are fast and reliable 



33 



34|. 



A simple model but which nevertheless presents all the essential behavior of bursting 



neurons is the Rulkov map 

mm 

Xn+l = ^ o + Vn, (33) 
Vn+l =yn- (y^n " (34) 

where x„ is the fast and ?/„ is the slow dynamical variable. The first variable has a dynamical 
behavior emulating the spiking-bursting activity of a neuron, depending on the parameter a, 
whereas the latter variable undergoes a slow evolution because of the small values taken on 
by the parameters a and /3, which model the action of external dc bias current and synaptic 



inputs on a given isolated neuron 201]. 



We suppose a transient chaotic behavior for the characteristic spiking of the fast variable 
x„, what is accomplished by choosing the values of the parameter a within the interval 
[4.1, 4.3] [Fig. [3]^a)]. The bursting timescale, on the other hand, comes about by the influence 
of the slow variable This can be understood by using a simple argument: since, from Eq. 
( 15^ . yn represents a small input on the fast variable dynamics its effect can be approximated 
by a constant value The resulting one- dimensional map, 

^-+1 = TT^ + ^^^^ 

can have either one, two, or three fixed points x\ 2 3, depending on the value of the input '0. 

As the latter approaches a critical value '^sn the fixed points x\ 2 (one stable and another 
unstable) undergo a saddle- node bifurcation, such that, for > however, the fixed 
points x\2 disappear [Fig. |2]. For values of t/^ > i\}c there is also a chaotic attractor that, 
provided tl}c < ip < ipsN, coexists with the stable fixed point attractor. Actually, at = ipc 
the chaotic attractor collides with the unstable fixed point xl and is destroyed through a 
boundary crisis [Fig. [2] (2lj| . The bursting regime then results from a hysteresis between the 
stable fixed point (quiescent evolution) and the chaotic oscillations (fast sequence of spikes). 

We consider a given burst to begin when the slow variable ?/„, which presents nearly 
regular saw-teeth oscillations, has a local maximum, in well-defined instants of time we call 
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■^-3 Vsn -2.5 -2 



FIG. 2. Bifurcation diagram for the map (j35p when a = 4.1. 2 3 fixed points, and we 

indicated the location of the saddle-node bifurcation (tpsN) and interior crisis (tpc)- 

Uk [Fig. [3]^b)]. We can define a phase describing the time evolution within each burst and 
varying from to 2tt as n evolves from to rifc+i: 

^{n) = 27ik + 2n ^ ~ . (36) 

Uk+l - Uk 

The duration of the chaotic burst, n^+i — Uk, depends on the variable x„ and fluctuates in 
an irregular fashion when x„ undergoes a chaotic evolution. There follows that the bursting 
phase rate also varies with time, such that we look at the bursting frequency defined by 

n = ^^M^^. (37) 

n 

B. Network structure 

When the Euclidean distance between neurons does not play a significant role, the cor- 
responding networks may be treated from the graph-theoretical point of view. However, 
once we regard those neurons as embedded in a three-dimensional lattice (the brain, where 
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FIG. 3. Time evolution of the (a) fast and (b) slow variables in the Rulkov map (|33p -(|34p for 
a = 4.1, a = P = 0.001. We also indicate by the local maxima of the slow variable, used to 
obtain a phase for the bursting dynamics. 

they are connected by axons and dendrites), it is convenient to use a lattice embedded in a 
Euclidean space jsS^ . 

Such higher-dimensional lattices can be very difficult to work with in terms of computer 
simulations, specially if long-range interactions are present and the number of neurons is 
large. However, good insights are expected to come from simpler models, which can never- 
theless retain some of the general characteristics of higher- dimensional lattices. It is from 
this point of view that we use one-dimensional lattices with neurons, each of them being 
described by the Rulkov map ( l33ll34l) . 

The main assumption we make is that the neurons are coupled through the release and 
absorption of a neurotransmitter, which is able to diffuse quickly in the medium in which 
the neurons are embedded. The release of the neurotransmitter is influenced by the fast 
dynamics, that is, once the action potential is spiking there are liberated neurotransmitter 
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molecules in the cell environment. These molecules are then absorbed by other neurons and 
are supposed to affect their fast dynamics, making them to spike in a slightly different way. 

This would be certainly a rather crude model to describe information transmission 
through neuron cells, for it does not incorporate other features of neuron electrical and 
chemical coupling. However, this model, in spite of its simplicity, suffices to exemplify in 
which sense does a coupling mediated by a diffusing substance can affect coherent rhythms 
and other collective phenomena, what is potentially interesting for other applications, like 
chemotaxis. 

The model of nonlocal coupling mediated by a diffusing substance, outlined in the previ- 
ous Section, can be straightforwardly translated to the language of coupled map lattices, for 
which both space and time are discrete. We will use the notation 'X.n^ for the state variable 
vector at the site j = 1, 2, . . . (in a one- dimensional lattice with periodic boundary condi- 
tions) and time = 0, 1, 2, . . .. The vector field F specifies now the map equations, and the 
interpretation of the remaining variables is the same as for oscillator chains. 

The equation for linear coupling of maps becomes 



the future coupled map lattice being obtained by replacing Xn by F(X„'^) in the coupling 
term. The normalization constant is still given by Eq. ( I25l) . 

Let each neuron to undergo a discrete time evolution governed by the Rulkov map (|33|) - 
(1341). for which M = 2 and 



where a, a and /3 are model parameters. Since there is always some biological diversity in a 
neuron assembly it is reasonable to choose slightly different values for the map parameters. 

In order to avoid chance correlations of bursting for uncoupled neurons we choose ran- 
domly the a parameter (influencing the fast dynamics) within the interval [4.1,4.3] with a 
uniform probability. The slow dynamics parameters a and /3 were all set up to the same 
(small) values, namely 0.001, since their possible differences are not likely to affect the results 
as the parameter a does. 




(38) 




(39) 
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Since we assume that the release and absorption of neurotransmitters affect the fast 
dynamics, we couple the Rulkov maps through the x-variable. For the sake of simplicity, we 
use a linear coupling. In some other applications, however, a future coupling would be better 
so as to keep the coupled variables within prescribed ranges, a feature that is warranted in 
linear couplings only if the coupling strength is weak enough. These assumptions amount 
to choose the following connection matrix 



A= , (40) 




where e is the coupling strength. The resulting linearly coupled Rulkov map lattice becomes 

N' 

- + y\i^ +eCY^ e--^^ \x^-^ + x^^'^] , (41) 

1=1 

yiili = y^^ - - /S, (42) 



1 + (x^f'V 



where we include a superscript to the a variables to denote the different values they may 
take throughout the network. 



V. BURSTING SYNCHRONIZATION OF COUPLED RULKOV NEURONS 

The coupled Rulkov map lattice fHT!) - fH2l) cannot exhibit a completely synchronized state, 

^(1) ^ ^(2) ^ . . . ^ AN) „,(1) ^ „,(2) ^ ... ^ .,{N) 

•^n •^n •^n ' Un tin tin 

For this to occur it would be necessary that the completely synchronized state be a pos- 
sible solution of Eqs. fHT]) - fH2]) . stable under infinitesimal perturbations along directions 
transversal to this state. One of the conditions for this is that the individual maps must 
be identical. This is a condition too stringent to be obeyed by realistic neuron models, for 
which there must be some degree of diversity in parameters. To mirror this fact, in our 
numerical simulations the values of the parameter a of each neuron were randomly chosen 
inside a given interval, provided we always have bursting. However, as a consequence, a 
completely synchronized state is not possible for our model. 

This does not necessarily mean that the neural network cannot present some degree of 
coherent behavior. The neuron bursting phases, for example, can synchronize through the 

15 



interaction provided by the coupling. We measure this effect by computing the mean field. 



N 
i=i 

If the neurons are weakly coupled, they burst at different times in a non-coherent fashion, 
and the mean field fluctuates irregularly with small amplitudes. Oppositely, if the neurons 
burst synchronously (i.e. they start bursting at approximately the same times) a nonzero 
mean field is formed and M„ presents regular oscillations of comparatively large amplitude. 
Only the slow dynamics becomes coherent as the neurons burst synchronously. The fast 
(spiking) dynamics remains incoherent and do not contribute to the mean field dynamics. 



which is kept close to a periodic regime 22 1 



When the neurons are more globally coupled (small 7) the mean field exhibits large- 
amplitude oscillations indicating that the neurons are bursting at approximately the same 
times [Fig. Hl^a)]. This can also be shown by inspecting the individual time evolution of the 
fast variable [Figs. Ilt^d) and (g)] where we select two neurons with different values of a, 
which start bursting at nearly the same time instants. Since such behavior is likely to be 
observed for all other neurons this corresponds to a coherent output for the entire network. 

On the other hand, when the neurons are more locally coupled, what is achieved by 
choosing larger values of 7, only the nearest neighbors contribute appreciably to the coupling, 
and the mean field displays only small-amplitude fluctuations, as illustrated by [Fig. 
Hl^c)]. In this case different neurons burst out of phase [Figs. ID^f) and (i)]. An intermediate 
value of 7 parameter shows a nearly coherent neuron bursting [Fig. H] (b)], for the neurons 
burst at slightly different times [Figs. Ht^e) and (h)]. This suggests that, as the coupling 
becomes increasingly global, there is a transition to bursting synchronization. 

In order to characterize bursting synchronization, however, the mean field is not precise 



enough, such that we use instead Kuramoto's order parameter 30 1 



1 ^ 

Zn = Rnexp (i<l>„) = — ^exp(i(/?(f)), (44) 

i=i 

where i?„ and $„ are the amplitude and angle, respectively, of the order parameter. If the 
neurons burst at exactly the same times, their phases superimpose coherently such that 
Rn = 1. By way of contrast, if the neurons burst in completely uncorrelated times the 
corresponding phases (pn^ would add to a near-zero value. Hence Rn can be used as a 
numerical diagnostic of bursting synchronization. 
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FIG. 4. Time evolution of the mean field (a) and the fast variables of two selected neurons [(d) and 
(g)] belonging to a linearly coupled lattice of = 251 Rulkov neurons with 7 = 0.005. (b,e,h) are 
the corresponding results for 7 = 0.0125, and (c,f,i) for 7 = 0.05. The remaining parameters are 
a = /3 = 0.001 and the parameter a is randomly chosen within the interval [4.1, 4.3], with coupling 
strength e = 0.1 

In Fig. |5]we plot the time-averaged order parameter magnitude < R > (after the tran- 
sients have died out) as a function of the coupling strength e, for different values of 7. For 
very small values of e the order parameters take on very small, yet nonzero, values. Hence 
the weakly coupled neurons are not likely to produce coherent bursting activity. The order 
parameter magnitude would not be equal to zero, even for very weak coupling, since the 
lattice size N is finite. As we consider larger lattices this value is expected to approach zero. 

As a general trend we see that a transition to synchronized bursting as the coupling 
strength increases. Specially for very small values of 7, this transition is similar to that 
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FIG. 5. Time-averaged order parameter as a function of the coupling strength for 7 equal to 0.025 
(circles), 0.020 (squares), 0.0125 (diamonds), 0.005 (triangles) and 0.002 (crosses). We remaining 
parameters are the same as in Fig. HI 



exhibited by the well-known the Kuramoto model (or globally coupled oscillators) [see Eq. 
(fT6|) ] [30]. In fact, the latter would correspond to the 7 — )■ limit of our model. Provided it 
is kept small enough, as 7 increases, this transition occurs for higher values of the coupling 
strength. However, for a further increase, this transition no longer seem to occur. 

These observations can also be drawn from analyzing Fig. [6|, where the time-averaged 
order parameter magnitude is plotted against the range parameter 7, for different coupling 
strengths. If the range parameter is large enough, no coupling strengths seem to produce 
any synchronization at all. On the other hand, if the range parameter is zero, bursting syn- 
chronization only occurs after a critical coupling strength, what is confirmed from inspecting 
the curves in Fig. [6] for 7 close to zero. 
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FIG. 6. Time-averaged order parameter as a function of the range parameter for e equal to 
0.0125 (circles), 0.0250 (squares), 0.5000 (diamonds), 0.1000 (triangles) and 0.2000 (crosses). We 
remaining parameters are the same as in Fig. HI 



VI. CONTROLLING BURSTING SYNCHRONIZATION 

As stated in the Introduction, bursting synchronization, in the context of biological neu- 
ron assemblies, can lead to undesirable and even pathological rhythms. Hence one may 
think of a control strategy to suppress bursting synchronization. Such a control procedure, 



although in a tria^ 
brain stimulation 



-and-error basis, has been implemented in neurosurgery through deep 



micro-electrodes are implanted in deep brain 
regions of a patient, like the subthalamic nucleus or globus pallidus, and a high-frequency 
(in the 100-120 Hz range) low- amplitude signal is applied 



We will consider two strategies for controlling bursting synchronization. One of them 
consist on the injection of a time-periodic electric signal of low amplitude (so as not to 
damage the neurons) and whose ultimate goal is the suppression of bursting synchronization, 
leading to normal neural rhythms. The second technique is the use of a time-delayed feedback 
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signal, with the same purpose. Each control method has its advantages in terms of the 
efficiency of the suppression of synchronized bursting. 

The effectiveness of each control procedure on reducing or suppressing synchronization 



can be measured by the suppression coefficient 



25| 



^=^/Va^' 

where X and Xf are the values of the mean field in the absence and presence of the control, 
respectively. A control scheme is ideally efficient when the variance of the controlled mean 
field vanish, irrespectively of its value without control, corresponding thus to an infinite 
value of S. As a general rule, the larger is the value of S, the more efficient will be the 
feedback on suppressing synchronization. 



A. Injection of a time-periodic signal 

We can adapt the map-based neural network model fHT]) - fH2|) so as to model the injection 
of a time-periodic electric signal of low amplitude. We have implemented this procedure by 
adding an external time-periodic signal to just one neuron j = S (the remaining neurons 
remaining unchanged). 

Moreover, in the same way we have proceeded in the previous Section, we have chosen 
slightly different values for the neuron parameter 4.1 < a*^*-* < 4.3. The slow timescale 
parameters a and (3 have been kept fixed for all neurons, since small variations on them 
have caused no noticeable effect. The controlled coupled map lattice is thus 



(j) 



N' 



1 + (xlfO e=i 
y^l = y\i^ - ax^ - /3, (47) 

where d and uj are the external signal amplitude and frequency, respectively, and 5j^s is the 
Kronecker delta. The neuron j = S, on which the control signal is applied, can be randomly 
chosen. The coupled equations are also straightforwardly modified on considering situations 
where more than one neuron is acted upon, with a signal with the same amplitude and 
frequency. 

We remark that the time-periodic signal is applied only at the variable representing the 
fast (x) neuron dynamics, since the slow {y) variable modulates the chaotic activity of x. 
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For this to occur even when the neurons are coupled, it is necessary that the external inputs 
on the slow variables be comparatively small, a condition that cannot be warranted if the 
coupling is made to occur also in both variables. 

1. Suppression efficiency 



possibility of bursting phase synchronization, as we dealt with in the previous Section, we 
shall also be concerned with bursting frequency synchronization, which is a slightly weaker 
phenomenon since it demands that only the time-rate be equal for a number of neurons, 
even when their phases may differ. We also consider frequency synchronization, in the con- 
text of the present model, as an undesirable feature, and we use it so as to investigate the 
controllability of the time-periodic signal with a given "external" frequency u {25! . 




In Fig. [Tl^a) we plot (in a colorscale) the suppression coefficient given by Eq. ( H5|) as a 
function of the amplitude d and frequency a; of a time-periodic control signal applied to the 
site S =1 of a network of = 111 Rulkov neurons with 4.1 < a^^^ < 4.3, a = /3 = 0.001, e = 
0.1, and 7 = 0.005 (near zero, hence in the most favorable case since the coupling is almost 
global). The coupled Rulkov map lattice presents a robust synchronization of bursting, since 
there is little effect of the external signal on suppressing synchronized bursting. 

It is possible to understand this robustness of bursting by regarding the synchronized 
coupled map lattice as a single oscillator with many degrees of freedom being perturbed by 
a time-periodic forcing. Hence we expect to see a mode-locking type of phenomenon, with 
many "Arnold-like" tongues representing values of amplitude for which the system locks in 
the frequency of the external signal. In Fig. [8] we show some Arnold tongues corresponding 
to four different values of the 7 parameter (represented by different colors). Since these 
tongues are of finite width, small changes in the signal amplitude cannot drive the system 
out of these tongues (unless we are close enough to the tongue border). 

A comparison of interest consists on keeping the signal frequency constant and plotting 
the suppression efficiency as a function of the signal amplitude d and the coupling strength 
e, for different values of the 7 parameter, as before [Fig. El^b)]. Increasing the signal 
amplitude does not lead to an appreciable enhancement of the efficiency on suppressing 




^^\0))/n stands for the corresponding frequency. Besides the 
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FIG. 7. (color online) Suppression coefficient as a function of (a) amplitude and frequency and (b) 
coupling constant and amplitude, for time-periodic control signal of a network of = 111 Rulkov 
neurons with 4.1 < a^^^ < 4.3, a = /3 = 0.001, 7 = 0.005. In (a) we used e = 0.1 and, in (b) we 
used uj = 0.1. The driving signal is applied to the site 5 = 1. 

bursting synchronization. Although not shown in Fig. Wljo), we did numerical simulations 
for other values of 7, confirming that the synchronized behavior is rather insensitive to 
changes in the signal amplitude. 



2. Bursting frequency locking 

Many of the features we have seen for a time-periodic signal can be explained by in- 
terpreting its action on the network as a frequency- locking phenomenon. Hence we shall 
describe some aspects of this influence that are similar to the mode-locking phenomena 



39| 



usually described in nonlinear oscillations perturbed by harmonic forcing 

We start by considering values of the coupling strength e for which the unperturbed 
lattice {d = 0) exhibits bursting synchronization, with frequencies il^^^ distributed in a 
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FIG. 8. (color online) Frequency-locking "Arnold tongues" for bursting frequency in the external 
signal parameter plane of amplitude vs. frequency, for 7 = 0.015 (black); 0.025 (red); 0.035 (blue), 
and 0.05 (brown). The remaining parameters are the same as in the previous figure. 

very narrow interval due to the distinct parameter values of each neuron. When the time- 
periodic signal is applied, however, the bursting frequencies ^l^^^ exhibit a locking with the 
external frequency u, what is represented in Fig. [9] as a horizontal plateau for the frequency 
mismatch Q^^'' — w, which is plotted against uj for all neurons belonging to the network. The 
concentration of these points reflect the existence of bursting synchronization, such that the 
latter is more suppressed as the thicker is the corresponding curve. 

In this sense the ability of suppress bursting synchronization is better for higher values 
of 7 [Fig. Mjo)] than for lower ones [Fig. [U^a)], since the curves in (b) are thicker than the 
curves in (a). The horizontal plateaus representing frequency locking increase in both cases 
with the signal amplitude d, which is a general property of mode locking. On inspecting 
Fig. [8] we see that, for different values of 7 parameter, the tongues arise from different values 
of the external signal frequency u. This "zero-amplitude" frequency Qo increases with the 
parameter 7, i.e. with couplings more and more of a locally than a globally nature, as shown 
by Fig. M 

However these tongues are different from those observed in simpler systems like the 
sine-circle map, for example: their widths increase monotonically from both sides but with 
different inclinations. Moreover, after some value of d, the increase seems to saturate where 
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FIG. 9. (color online) Frequency mismatch of bursting neurons versus the external driving fre- 
quency for a lattice with = 51 neurons, 4.1 < a^^^ < 4.3, a = j3 = 0.001, e = 0.1, and (a) 
7 = 0.025; (b) 7 = 0.050. In both cases, the driving signal applied at the site 5=1 with ampli- 
tude d = (black points), 0.05 (red points), and 0.15 (blue points). These points are plotted for 
all network sites. 

the inclination was higher, whereas it continues to increase at the same rate where the 
inclination was lower. A partial explanation of this effect is that the frequency- locking 
interval is different with respect to the point where c? = 0. 

On inspecting again Fig. |9] we see that the zero-amplitude frequency VLq of each interval 
is the intersection of the line black points (for which c? = 0) and the horizontal line (where 
= cj). The interval is thus wider at the righthandside of ojq than at its lefthandside. 
Hence it is useful to define not only the total interval width Aw but also its lefthand length 
5uj. If the tongues were nearly symmetrical, we would have 5uj ~ /S.VL/2. 

In Fig. [11] we plot both widths as a function of the signal amplitudes for different values 
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FIG. 10. Zero-amplitude frequency (from where Arnold tongues arise in the parameter plane) as 
a function of the 7 parameter for the same parameters as in the previous figures. 



of 7. The overall behavior of them is similar, but with some noteworthy differences: for very 
small signal amplitudes the tongue widths are nearly independent of the 7 value, signaling 
that the influence here is more from the individual dynamics than from the coupling itself. 
The latter manifests only after some threshold, and the widths are generally smaller for 
higher 7-values than for smaller ones. 

The asymmetric character of the tongues also depends on the coupling parameter 7. Let 
us take a fixed 7, such as 0.0125 (represented as circles in Fig. [H]). As the signal amplitude 
increases both the total width Au [Fig. [TTT a)] as the righthand width 6u [Fig. [TTT b)] increase 
very mildly with d. However, while the former width takes on values around 5.5 x 10~^, 
for large d, the latter is restricted to ~ 4.5 x 10~^, which is nonetheless an asymmetric 
situation. On the other hand, if we consider a larger 7, such as 0.05 (diamonds in Fig. [TTjl . 
the total width Au is of the order of 4 x 10~^ [Fig. ITlT a)]. whereas the righthand width 
6u is nearly half of this value [Fig. [TTT b)]. Hence for larger 7 we have more sjTiimetrical 
Arnold tongues than for smaller 7, suggesting that the asymmetry of the Arnold tongues is 
chiefly a coupling effect, at least as the effective range is concerned. 

A qualitative explanation for the asymmetry of the mode-locking tongues has been pro- 
vided by Ivanchenko et al. 22|: an imposed positive signal precipitates a burst into a 
quiescent regime, and delays it when the signal is negative. If the driving frequency is 
higher than that of the mutually synchronized network the periodic signal will fasten the 
oscillations of the neurons. On the other hand, when the driven neuron starts a burst, there 
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FIG. 11. Frequency locking tongue widths as a function of the signal amplitude for (a) total width 
Ao; and (b) righthand width 6uj (with respect to the zero-amplitude frequency VLq). In both panels 
we have 7 = 0.0125 (circles); 0.035 (squares); and 0.050 (diamonds). The remaining parameters 
are the same as in the previous figures. 



is an abrupt change in the mean field perceived by all neurons in the lattice, pushing them to 
a quiescent state. As a result, higher frequencies would give better synchronization effects. 



B. Delayed feedback control 



The suppression of bursting synchronization through a time-delayed feedback signal has 
been proposed by Rosenblum and Pikowsky [2J, l25|]. The lattice coupling is represented by 
a term eX„, which includes the mean field X„ given by Eq. (143|) . According to Ref. 2J] we 
consider two procedures of feedback control, with respect to their dependences on the mean 
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field, described by 



(0 








1 + (x« 


Vn+l 





+ i/« + eXn + EfX^., - e'Xr,, (48) 
(z = l,2,...iV) (49) 
where 

N' 

X„ = C ^ e-^^^ + x(f +^)] (50) 

1=1 

and such that we have: (i) direct feedback, which takes into account the current mean field 
and its value r iterations before, X„_t-, with coupling intensity £f, in such a way that e' = 0; 
(ii) differential feedback, for which the controlling term is the difference between the current 
and the time-delayed mean fields, with e' = sj. 



1 . Differential feedback 

An example of differential feedback control to suppress bursting synchronization is given 
in Fig. [121 where we plot the suppression coefficient as a function of the control parameters 
Ef and r. The latter are normalized according to the coupling strength e = 0.1 and the 
mean bursting period T = 200, respectively. 

In contrast with the case of a time-periodic signal, now we observe the formation of 
suppression domains, or regions with large (or small) values of S, in the control parameter 
plane. As a general trend (i.e., for all 7 values considered), we observe that the domains 
of high (low) suppression occur for positive (negative) values of the control strength, these 
values being symmetric for a given value of the time delay r. 

The domains are roughly centered at values of r which are integer multiples of the bursting 
period T. Moreover, the suppression domains (both high and low) are periodic in the time 



delay. Both features 
with global coupling 



lave been already observed by Rosenblum and Pikowsky in networks 
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25l |. In our case, we observe an infiuence of the parameter 7: as it 
increases (meaning that our coupling becomes more local) the suppression domains become 
centered in progressively smaller values of r and the period in r also decreases. On the other 
hand, the maximum efficiency seems not to depend on the 7, since it varies between 1.04 
[Fig. [I2]^a)] to 1.10 [Fig. [I2]^c)], a mere ~ 5% of increase. 
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FIG. 12. (color online) Suppression coefficient as a function of the normalized control parameters 
(strength and time delay) for differential feedback control of a network of = 111 Rulkov neurons 
with r = 200, e = 0.1 and (a) 7 = 0.005; (b) 0.0125; (c) 0.025. The remaining parameters are the 
same as in the previous figures. 

2. Direct feedback 

Now we consider in Fig. [13] an example of direct feedback control, with S plotted against 
normalized ej and r, as before. We still observe suppression domains, but with a marked 
contrast with the differential feedback case. Firstly the domains of low suppression are 
clearly more often observed than those of high suppression, showing a striking asymmetry. 
The domains of low suppression present three noteworthy features: (i) they occur both for 
positive and negative values of the control strength; (ii) they are periodic in r; (iii) they 
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FIG. 13. (color online) Suppression coefficient as a function of the normalized control parameters 
(strength and time delay) for direct feedback control of a network of Rulkov neurons with (a) 
7 = 0.005; (b) 0.0125; (c) 0.025. The remaining parameters are the same as in the previous figure. 

are roughly centered at multiple integers of r/T; and (iv) they are dimerized, i.e. a domain 
of low suppression for positive is followed by another one for negative Ef, after circa one 
period T. 

The location of the low efficiency domains seems not to be noticeably affected by 7, but we 
do observe such a domain for very small r only in the case of globally coupled lattice (smaller 
7). A curious observation is that there is a tiny domain of very high efficiency (higher than 
any other domain in the differential control) for negative Ef, this feature observed for all 7 
values considered. 

This does not mean, however, that the direct feedback is worse than the differential 
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feedback, as far as its suppression efficiency is concerned. Comparing Figs. [12] and [13] we 
observe that the regions of moderate efficiency of the latter (between the low efficiency 
domains) have values of S comparable with the domains of high efficiency of the differential 
control. However, as a whole there are more regions of high efficiency in the direct control 
than in differential one, the difference being that the high efficiency regions are spread out 
over the parameter plane in the direct feedback case. 

VII. CONCLUSIONS 

In many biological contexts the interaction between dynamically active cells is mediated 
by a chemical diffusing through the inter-cell medium. This is the case, for example, of 
neuron synaptical connections, which are accomplished through neurotransmitters. There 
are two timescales involved in this kind of situations, since the timescale of the cell oscillation 
is generally different from the timescale related to the diffusion process. Hence a complete 
treatment of this kind of coupling would demand the simultaneous solution of the diffusion 
equation and the system of coupled oscillator equations. 

However, if the timescale related to the diffusion is much smaller than the characteristic 
oscillator periods we can make an adiabatic approximation and write down a closed-form 
expression for the nonlinearly coupled oscillator chain. For the one-dimensional case it results 
in a coupling term whose strength decreases with the lattice distance in an exponential 
fashion. An advantage of this procedure is that one can describe the commonly used global 
(all-to-all) and local (nearest-neighbor) limiting forms of the nonlocal coupling, by 

varying the diffusion length that characterizes the exponential decay of interaction. 

In this work we explored some of these features by using a simplified model for neuronal 
network consisting of map-based units displaying spiking-and-bursting behavior. Although 
the spiking (fast) dynamics is chaotic, the bursting (slow) dynamics presents coherent fea- 
tures due to the neuron interaction, one of them being their capability of synchronize the 
beginning of the bursting regime. 

The transition to bursting synchronization was found to depend on the nonlocal coupling 
features. In particular, as the inverse coupling length 7 is decreased it becomes more dif- 
ficult to obtain synchronization, for a fixed coupling strength. This is compatible with the 
observation that locally coupled neurons are less amenable to exhibit such collective effects 
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like bursting synchronization than globally coupled ones. 

This is particularly interesting from the point of view of potential disorders (like abnormal 
rhythms associated with epilepsy and Parkinson disease) related to them. In order to control 
those undesirable bursting synchronized rhythms we have used two strategies for suppressing 
bursting synchronization: a time-periodic signal and a delayed feedback one. In the former 
procedure a time-periodic harmonic and low-intensity signal injected on a randomly chosen 
neuron. This external signal is capable to change the bursting frequencies of interacting 
neurons so as to drive them out of the synchronized behavior. The inclusion of a time- 
periodic signal causes the synchronized neurons to lock into an Arnold-like tongue with a 
well-defined width. 

The common bursting frequency was found to increase with the inverse coupling length. 
Moreover, the frequency locking tongues were found to be asymmetric with respect to this 
common locking frequency. As a consequence, if one injects an external signal with frequency 
higher than the upper limit of the frequency-locking tongue we can obtain desynchronization, 
as required. Both the amplitude and the asymmetrical features of the tongues were found 
to depend on the nonlocal character of the interaction between neurons. 

In terms of the efficiency of the bursting synchronization suppression, we did not find an 
appreciable effect of an external time-periodic signal of constant amplitude and frequency, 
even after wide changes in these parameters. Accordingly, we sought for a non-constant 
scheme of control, preferably with an amplitude tailored to suit the actual needs of sup- 
pressing synchronized behavior. This was accomplished by the injection of a time-delayed 
feedback signal. In contrast with the previous control scheme, for this feedback signal we 
observed formation of suppression domains for time delays roughly centered at multiples of 
the bursting period, with a well-defined periodicity. 

Two feedback control types were used: a differential scheme, with a relative symmetry 
between domains of low and high suppression, and a direct scheme, with an apparent bias 
towards low suppression domains. On the other hand, we found that the direct scheme 
seems superior to both a time-periodic signal and a differential feedback signal in the sense 
that the suppression efficiency is higher over the control parameter plane, even though we 
do not observe many domains of control for this case. Hence, from a "practical" point of 
view, a direct feedback delayed signal would yield better results in terms of the suppression 
of undesired (and pathological) synchronized rhythms. 
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The extension of the present formahsm to higher-dimensional lattices is quite straightfor- 
ward, but the coupling kernel is somewhat different and the summation through all neigh- 
borhoods may become a computationally difficult task. However, the general features of 
such a kind of coupling are already present in the simpler one- dimensional case. 
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